module parameter_setting

    implicit none

    ! model parameters

    ! parameters given by a matlab code ###############################################################3

    double precision :: parameters(13)
    double precision :: parameters_ss(16)
    double precision :: parameters_integer(8)
    double precision :: parameters_solution(2)
    double precision :: init_guess(6)
    integer :: init_type
    
    double precision :: policy_parameters(3), tau_l, thre_l, thre_n
    double precision :: welfare_states(3), n_stss, k_stss, rb_stss
    double precision :: welfare_integer(2)
    integer :: length, repeat
    integer, allocatable :: pol_flag(:)

    double precision :: beta
    double precision :: nu
    double precision :: alpha
    double precision :: delta
    double precision :: lambda
    double precision :: chi0
    double precision :: chi1
    double precision :: rho_a
    double precision :: sigma_a
    double precision :: nd

    double precision :: psi
    double precision :: n0
    double precision :: gamma
    double precision :: sigma_rk
    
    double precision :: c_ss,y_ss,i_ss,k_ss,h_ss,n_ss,l_ss,rb_ss,rkhats_ss,a_ss,xi_ss,pr_ss

    ! simulation
    double precision, allocatable :: av(:),kv(:),nv(:),rbv(:),xiv(:)
    double precision, allocatable :: cv(:),yv(:),iv(:),hv(:),lv(:),rhseev(:),pibv(:)
    double precision :: rkhat,rkhats,rkhat_next,rkhats_next
    double precision :: erk_next,pibv_next
    double precision :: xi_next,c_next,h_next,y_next,rhsee_next,rv_next,a_next
    double precision :: n_next,i_next,k_next,l_next
    double precision :: z_temp
    double precision :: z_temp_log, cdf_log
    
    double precision, allocatable :: x0(:,:),x(:),x0_next(:,:),x_next(:),x0_reg(:,:),x_reg(:)
    double precision, allocatable :: rhsee_each(:),rhsee_derived(:)
    double precision :: dummy1,dummy2,rhsee_temp,cdf,pdf
    double precision :: mean_rhsee_norun,mean_rhsee_run,mean_rhsee_nega
    double precision :: mean_rbar_norun ,mean_rbar_run ,mean_rbar_nega, mean_rbar_pol
    double precision, allocatable :: rbar_derived(:)
    double precision :: rbar_derived_temp

    ! polynomials
    integer :: nstv,degp,npol ! number of state variables, degree of polynomial, size of polynomial matrix
    double precision, allocatable :: expmx(:,:),expmx_column(:),state(:)
    
    double precision, allocatable :: eta1_rhsee(:),eta2_rhsee(:),eta3_rhsee(:),eta4_rhsee(:)
    double precision, allocatable :: eta1_rhsee_new(:,:),eta1_rhsee_old(:,:)
    double precision, allocatable :: eta2_rhsee_new(:,:),eta2_rhsee_old(:,:)
    double precision, allocatable :: eta3_rhsee_new(:,:),eta3_rhsee_old(:,:)
    double precision, allocatable :: eta4_rhsee_new(:,:),eta4_rhsee_old(:,:)
    
    double precision, allocatable :: eta1_rbar(:),eta2_rbar(:),eta3_rbar(:),eta4_rbar(:)
    double precision, allocatable :: eta1_rbar_new(:,:),eta1_rbar_old(:,:)
    double precision, allocatable :: eta2_rbar_new(:,:),eta2_rbar_old(:,:)
    double precision, allocatable :: eta3_rbar_new(:,:),eta3_rbar_old(:,:)
    double precision, allocatable :: eta4_rbar_new(:,:),eta4_rbar_old(:,:)

    integer :: con_max(1)
    double precision :: con_max_new(1,1), con_max_old(1,1)
    
    ! trapezoid integration
    integer :: n_trapz
    double precision, allocatable :: z_trapz(:)

    ! policy function iteration and convergence
    integer :: npf
    double precision :: tol, adj
    double precision, allocatable :: con_lev(:), rsq(:), shocks(:)
    integer :: t_max, iter_max, t_drop
    integer :: t, i, j, f, g, z, count

    ! values for non-linear solver
    integer :: nnls
    double precision, parameter :: tol_fun = 1d-8   
    double precision :: guess(1),fval(1)
!    double precision, allocatable :: guess(:),fval(:)
    integer :: info
    integer :: nonsol

    ! regression to update polynomial coefficients
    integer :: num_run, num_norun, num_nega, num_pol
    double precision, allocatable :: x_run(:,:),x_norun(:,:),x_nega(:,:),x_pol(:,:)
    double precision, allocatable :: xx_norun(:,:),xx_inv_norun(:,:)
    double precision, allocatable :: xx_run  (:,:),xx_inv_run  (:,:)
    double precision, allocatable :: xx_nega (:,:),xx_inv_nega (:,:)
    double precision, allocatable :: xx_pol  (:,:),xx_inv_pol  (:,:)
    
    double precision, allocatable :: y_rhsee_run(:,:),y_rhsee_norun(:,:),y_rhsee_nega(:,:),y_rhsee_all(:)
    double precision, allocatable :: y_rhsee_pol(:,:)
    double precision, allocatable :: y_rbar_run(:,:) ,y_rbar_norun(:,:) ,y_rbar_nega(:,:) ,y_rbar_all(:)
    double precision, allocatable :: y_rbar_pol(:,:)
    double precision, allocatable :: xy_rhsee_norun(:,:),xy_rhsee_run(:,:),xy_rhsee_nega(:,:)
    double precision, allocatable :: xy_rhsee_pol(:,:)
    double precision, allocatable :: xy_rbar_norun(:,:) ,xy_rbar_run(:,:),xy_rbar_nega(:,:)
    double precision, allocatable :: xy_rbar_pol(:,:)
    
    double precision, allocatable :: y1_temp_run(:,:) ,y1_temp_norun(:,:) ,y1_temp_nega(:,:) ,y1_temp_pol(:,:)
    double precision, allocatable :: y2_temp_run(:,:) ,y2_temp_norun(:,:) ,y2_temp_nega(:,:) ,y2_temp_pol(:,:)
    double precision, allocatable :: x_temp_run(:,:)  ,x_temp_norun(:,:)  ,x_temp_nega(:,:)  ,x_temp_pol(:,:)

    ! to import data
    character(len=20) :: filename               

    double precision, allocatable :: solved_rbv(:)
    
    ! EE error
    double precision, allocatable :: cv_ee(:),error(:)
    double precision :: error_avg

end module

! main program ##########################################################################################   

program baseline_model

    ! modules
    use parameter_setting
    use my_toolbox

    implicit none
    
    external residual_fn

    ! import parameters from text files #################################################################

    open (10, file='from_matlab_to_fortran/parameters.txt', status='old', action='read')
       read(10,*) parameters
    close(10)
    beta     = parameters(4)
    nu       = parameters(5)
    alpha    = parameters(6)
    delta    = parameters(7)
    lambda   = parameters(8)
    chi0     = parameters(9)
    chi1     = parameters(10)
    rho_a    = parameters(11)
    sigma_a  = parameters(12)
    nd       = parameters(13)

    open (10, file='from_matlab_to_fortran/ss_values.txt', status='old', action='read')
       read(10,*) parameters_ss
    close(10)
    c_ss      = parameters_ss(1)
    y_ss      = parameters_ss(2)
    i_ss      = parameters_ss(3)
    k_ss      = parameters_ss(4)
    h_ss      = parameters_ss(5)
    n_ss      = parameters_ss(6)
    l_ss      = parameters_ss(7)
    rb_ss     = parameters_ss(8)
    rkhats_ss = parameters_ss(9)
    a_ss      = parameters_ss(10)
    xi_ss     = parameters_ss(11)
    pr_ss     = parameters_ss(12)
    psi       = parameters_ss(13)
    n0        = parameters_ss(14)
    gamma     = parameters_ss(15)
    sigma_rk  = parameters_ss(16)

    open (11, file='from_matlab_to_fortran/parameters_integer.txt', status='old', action='read')
        read(11,*) parameters_integer
    close(11)
    npf      = nint(parameters_integer(1)) ! number of policy functions to be approximated
    nstv     = nint(parameters_integer(2)) ! number of state variables
    degp     = nint(parameters_integer(3)) ! degree of polynomial
    npol     = nint(parameters_integer(4)) ! size of polynomial matrix
    iter_max = nint(parameters_integer(5))
    n_trapz  = nint(parameters_integer(6)) ! number of grids for integration approximation
    t_max    = nint(parameters_integer(7))
    t_drop   = nint(parameters_integer(8))

    open (11, file='from_matlab_to_fortran/parameters_solution.txt', status='old', action='read')
        read(11,*) parameters_solution
    close(11)
    tol = parameters_solution(1)
    adj = parameters_solution(2)

    allocate(expmx(nstv,npol))
    allocate(expmx_column(npol*nstv))
    open (11, file='from_matlab_to_fortran/polynomial_matrix.txt', status='old', action='read')
        read(11,*) expmx_column
    close(11)
    f=0
    do i=1,nstv
        do j=1,npol
            f=f+1
            expmx(i,j) = expmx_column(f)
        end do
    end do

    allocate(z_trapz(n_trapz))
    open (11, file='from_matlab_to_fortran/z_trapz.txt', status='old', action='read')
        read(11,*) z_trapz
    close(11)

    allocate(shocks(t_max+1))
    open (11, file='from_matlab_to_fortran/shocks.txt', status='old', action='read')
        read(11,*) shocks
    close(11)

    allocate(eta1_rhsee(npol),eta2_rhsee(npol),eta3_rhsee(npol),eta4_rhsee(npol))
    allocate(eta1_rhsee_new(npol,1),eta2_rhsee_new(npol,1),eta3_rhsee_new(npol,1),eta4_rhsee_new(npol,1))
    allocate(eta1_rhsee_old(npol,1),eta2_rhsee_old(npol,1),eta3_rhsee_old(npol,1),eta4_rhsee_old(npol,1))
    allocate(eta1_rbar(npol),eta2_rbar(npol),eta3_rbar(npol),eta4_rbar(npol))
    allocate(eta1_rbar_new(npol,1),eta2_rbar_new(npol,1),eta3_rbar_new(npol,1),eta4_rbar_new(npol,1))
    allocate(eta1_rbar_old(npol,1),eta2_rbar_old(npol,1),eta3_rbar_old(npol,1),eta4_rbar_old(npol,1))

    allocate(con_lev(npf),rsq(npf))

    open (11, file='from_fortran_to_matlab/solved_eta1_rhsee.txt', status='old', action='read')
    read(11,*) eta1_rhsee
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta2_rhsee.txt', status='old', action='read')
    read(11,*) eta2_rhsee
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta3_rhsee.txt', status='old', action='read')
    read(11,*) eta3_rhsee
    close(11)

    open (11, file='from_fortran_to_matlab/solved_eta1_rbar.txt', status='old', action='read')
    read(11,*) eta1_rbar
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta2_rbar.txt', status='old', action='read')
    read(11,*) eta2_rbar
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta3_rbar.txt', status='old', action='read')
    read(11,*) eta3_rbar
    close(11)
    open (11, file='from_fortran_to_matlab/solved_eta4_rbar.txt', status='old', action='read')
    read(11,*) eta4_rbar
    close(11)
    
    allocate(av(t_max+1),nv(t_max+1),kv(t_max+1),rbv(t_max+1),xiv(t_max+1))
    allocate(cv(t_max+1),yv(t_max+1),iv(t_max+1),lv(t_max+1),hv(t_max+1),rhseev(t_max+1),pibv(t_max+1))
    allocate(pol_flag(t_max+1))

    allocate(x_run(t_max,npol),x_norun(t_max,npol),x_nega(t_max,npol),x_pol(t_max,npol))
    allocate(xx_inv_run(npol,npol),xx_inv_norun(npol,npol),xx_inv_nega(npol,npol),xx_inv_pol(npol,npol))
    
    allocate(y_rhsee_run(t_max,1),y_rhsee_norun(t_max,1),y_rhsee_nega(t_max,1),y_rhsee_all(t_max),y_rhsee_pol(t_max,1))
    allocate(y_rbar_run(t_max,1) ,y_rbar_norun(t_max,1) ,y_rbar_nega(t_max,1) ,y_rbar_all(t_max) ,y_rbar_pol(t_max,1))
    allocate(xy_rhsee_run(npol,1),xy_rhsee_norun(npol,1),xy_rhsee_nega(npol,1),xy_rhsee_pol(npol,1))
    allocate(xy_rbar_run(npol,1) ,xy_rbar_norun(npol,1) ,xy_rbar_nega(npol,1) ,xy_rbar_pol(npol,1))

    allocate(y1_temp_run(t_max,1),y1_temp_norun(t_max,1),y1_temp_nega(t_max,1),y1_temp_pol(t_max,1))
    allocate(y2_temp_run(t_max,1),y2_temp_norun(t_max,1),y2_temp_nega(t_max,1),y2_temp_pol(t_max,1))
    allocate(x_temp_run(t_max,npol) ,x_temp_norun(t_max,npol) ,x_temp_nega(t_max,npol) ,x_temp_pol(t_max,npol))

    allocate(state(nstv))
    allocate(x0(nstv,npol),x(npol),x0_next(nstv,npol),x_next(npol),x0_reg(nstv,npol),x_reg(npol))
    
    allocate(rhsee_each(n_trapz))
    allocate(rhsee_derived(t_max))
    allocate(rbar_derived(t_max))

    allocate(solved_rbv(t_max))
    open (11, file='from_matlab_to_fortran/solved_rbv.txt', status='old', action='read')
        read(11,*) solved_rbv
    close(11)
    
    open (11, file='from_matlab_to_fortran/policy_parameters.txt', status='old', action='read')
        read(11,*) policy_parameters
    close(11)
    tau_l  = policy_parameters(1)
    thre_l = policy_parameters(2)
    thre_n = policy_parameters(3)

    open (11, file='from_matlab_to_fortran/welfare_states.txt', status='old', action='read')
        read(11,*) welfare_states
    close(11)
    n_stss  = welfare_states(1)
    k_stss  = welfare_states(2)
    rb_stss = welfare_states(3)

    open (11, file='from_matlab_to_fortran/welfare_integer.txt', status='old', action='read')
        read(11,*) welfare_integer
    close(11)
    length = nint(welfare_integer(1))
    repeat = nint(welfare_integer(2))

    allocate(cv_ee(t_max+1),error(t_max+1))

    ! iterate until policy function converges ########################################################

!    count = 1
!    do count = 1, iter_max

        ! set initial states
        av(1)  = a_ss
        nv(1)  = n_stss
        kv(1)  = k_stss
        rbv(1) = rb_stss
        lv(1)  = kv(1)/nv(1)
        
        pol_flag(:) = 0

        ! simulate for t_max periods
        do t = 2, t_max

        ! Variables determined independent of bank runs
        av(t)  = rho_a*av(t-1) + sigma_a*shocks(t-1)
        
        hv(t)  = ((1d0-alpha)/psi*exp(av(t))*kv(t-1)**(alpha))**(1d0/(alpha+1d0/nu))
        yv(t)  = exp(av(t)) * kv(t-1)**(alpha)*hv(t)**(1-alpha)
        
        rkhats = log ( rbv(t-1) * (1d0-1d0/lv(t-1)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
        rkhat  = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                + (1d0+nu)/(alpha*nu+1d0)*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t-1))
        xiv(t)  = exp(rkhat)/exp(rkhats)
        
        ! Up to this point, whether policy is binding or not is irrelevant.
        ! Whether cv should be obtained using rhseev depends on whether policy is binding or not.
        ! Proceed assuming policy is not binding, then later check.

        state(1) = log(kv(t-1))
        state(2) = log(rbv(t-1))
        state(3) = log(nv(t-1))
        state(4) = log(xiv(t))
        do i=1,nstv
            do j=1,npol
                x0(i,j) = state(i)**expmx(i,j)
            end do
        end do
        x  = x0(1,:)*x0(2,:)*x0(3,:)*x0(4,:)
        if ( xiv(t) >= 1d0 ) then
            pibv(t)  = ( exp(rkhat) + 1d0-delta )*lv(t-1)*nv(t-1) - rbv(t-1)*(lv(t-1)-1d0)*nv(t-1) - nv(t-1)
            if ( pibv(t) > 0d0 ) then
                nv(t) = chi1*( chi0*pibv(t) + nv(t-1) ) + n0
                rhseev(t) = exp( sum(x*eta1_rhsee) )
                rbv(t)    = exp( sum(x*eta1_rbar ) )
            else
                nv(t) = chi1*(      pibv(t) + nv(t-1) ) + n0
                rhseev(t) = exp( sum(x*eta3_rhsee) )
                rbv(t)    = exp( sum(x*eta3_rbar ) )
            end if
        else
            pibv(t)  = ( exp(rkhat) + 1d0-delta )*lv(t-1)*nv(t-1) &
                     - rbv(t-1)*(1d0+lambda*(1d0-gamma))*(lv(t-1)-1d0)*nv(t-1) - nv(t-1)
            nv(t) = n_ss*(1d0 - nd/100d0)
            rhseev(t) = exp( sum(x*eta2_rhsee) )
            rbv(t)    = exp( sum(x*eta2_rbar ) )
        end if
                
        cv(t) = rhseev(t)**(-1d0) + psi*hv(t)**(1d0+1d0/nu)/(1d0+1d0/nu)
        iv(t) = yv(t) - cv(t)
        kv(t) = (1d0-delta)*kv(t-1) + iv(t)
        lv(t) = kv(t)/nv(t)
        
        erk_next = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                 + (1d0+nu)/(alpha*nu+1d0)*rho_a*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t))
        
        ! To obtain deposit rate R_bar, solve bank's FOC w.r.t. levarage
        guess(1) = rbv(t) !log(erk_next) ! log(rbv(t-1))
        call hybrd1 ( residual_fn, 1, guess, fval, tol_fun, info, t )
        rbv(t)   = guess(1)
        
        if ( info /= 1 ) then
            nonsol = nonsol + 1
            write(*,*) 'Equation not solved at period',t
        end if

        rkhats_next = log ( rbv(t) * (1d0-1d0/lv(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
        z_temp     = (rkhats_next - erk_next) / sigma_rk
        z_temp_log = (erk_next + sigma_rk**2d0 - rkhats_next) / sigma_rk
        call normp( z_temp    , cdf    , dummy2, pdf    )    ! obtain cdf and pdf
        call normp( z_temp_log, cdf_log, dummy1, dummy2 )    ! obtain cdf
        rbar_derived(t) = exp(erk_next + sigma_rk**2d0/2d0) * cdf_log / (1d0-cdf) + (1d0-delta) &
                        - lambda * ( rbv(t)**2d0 * (1d0-gamma) * (1d0+lambda*(1d0-gamma)) * (lv(t)-1d0)/lv(t)**2d0 ) &
                                 / ( rbv(t) * (1d0-1d0/lv(t))  * (1d0+lambda*(1d0-gamma)) - 1d0+delta ) /sigma_rk * pdf/(1d0-cdf)
        
        ! cv,iv,kv,lv,erk_next,rbv need to be updated if policy is binding.
        ! Check whether policy is binding, and update cv,iv,kv,lv,erk_next,rbv if binding.
!        if ( xiv(t) >= 1d0 ) then
        if ( pibv(t) >= 0d0 ) then
        if ( lv(t) > thre_l .and. nv(t) > thre_n ) then
            pol_flag(t) = 1
            lv(t) = (1d0-tau_l)*lv(t)
            kv(t) = lv(t)*nv(t)
            iv(t) = kv(t) - (1d0-delta)*kv(t-1)
            cv(t) = yv(t) - iv(t)
            erk_next = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                     + (1d0+nu)/(alpha*nu+1d0)*rho_a*av(t) - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t))
            rbv(t)   = exp( sum(x*eta4_rbar) )
            ! rbv(t) is determined below so that Euler equation is satisfied.
        end if
        end if

        if ( pol_flag(t) == 0 ) then
        ! Numerical integration when policy is not imposed #############################################################
        ! Compute right-hand side of Euler equation using numerical integration ########################################
        ! Note that rhsee_derived(t=1:t_max) corresponds to E1(uc(2)),E2(uc(3)),... ####################################
            
        do z=1,n_trapz
        
            ! For each loop, take the same steps as above to compute endogenous variables.
            ! Each loop corresponds to a different xi shock.
            ! Use the variables above if xi shock does not affect their values.
            
            a_next  = rho_a*av(t) + sigma_a*z_trapz(z)
            h_next  = ((1d0-alpha)/psi*exp(a_next)*kv(t)**(alpha))**(1d0/(alpha+1d0/nu))
            y_next  = exp(a_next) * kv(t)**(alpha)*h_next**(1-alpha)
        
            rkhats_next = log ( rbv(t) * (1d0-1d0/lv(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
            rkhat_next  = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                        + (1d0+nu)/(alpha*nu+1d0)*a_next - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t))
            xi_next     = exp(rkhat_next)/exp(rkhats_next)
            
            state(1) = log(kv(t))
            state(2) = log(rbv(t))
            state(3) = log(nv(t))
            state(4) = log(xi_next)
            do i=1,nstv
                do j=1,npol
                    x0_next(i,j) = state(i)**expmx(i,j)
                end do
            end do
            x_next  = x0_next(1,:)*x0_next(2,:)*x0_next(3,:)*x0_next(4,:)
            
            if ( xi_next >= 1 ) then
                
                ! Derive l_next and n_next, then check if policy is binding.
                ! If binding, l_next=thre_l and c_next is determined accordingly.
                ! If not binding, use rhsee_next.
                
                pibv_next  = ( exp(rkhat_next) + 1d0-delta )*lv(t)*nv(t) - rbv(t)*(lv(t)-1d0)*nv(t) - nv(t)
                
                if ( pibv_next > 0d0 ) then ! if profit is positive
                    
                    rhsee_next = exp( sum(x_next*eta1_rhsee) )
                    c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
    !                rv_next    = rk_next * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rbv(t)
                    
                    n_next = chi1*( chi0*pibv_next + nv(t) ) + n0
                    i_next  = y_next - c_next
                    k_next  = (1d0-delta)*kv(t) + i_next
                    l_next  = k_next/n_next
                    
                    if ( l_next>thre_l .and. n_next>thre_n ) then
                        l_next = (1d0-tau_l)*l_next
                        k_next = l_next*n_next
                        i_next = k_next - (1d0-delta)*kv(t)
                        c_next = y_next - i_next
                        rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rbv(t)
                    end if
                
                else ! if profit is negative
                    
                    rhsee_next = exp( sum(x_next*eta3_rhsee) )
                    c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
    !                rv_next    = rk_next * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rbv(t)
                    
                    n_next = chi1*(      pibv_next + nv(t) ) + n0
                    i_next  = y_next - c_next
                    k_next  = (1d0-delta)*kv(t) + i_next
                    l_next  = k_next/n_next
                    
!                    if ( l_next>thre_l .and. n_next>thre_n ) then
!                        l_next = thre_l
!                        k_next = l_next*n_next
!                        i_next = k_next - (1d0-delta)*kv(t)
!                        c_next = y_next - i_next
!                        rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rbv(t)
!                    end if

                end if      
                
            else ! if run happens
                
                rhsee_next = exp( sum(x_next*eta2_rhsee) )
                c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
                rv_next    = ( exp(rkhat_next) + 1d0-delta ) * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rv_next
                
            end if

            ! multiply density, and cut off the both tails outside of z_trapz
            call normp( z_trapz(z), dummy1, dummy2, pdf )    ! obtain pdf at z_trapz(z)
            call normp( z_trapz(1), cdf, dummy1, dummy2 )    ! obtain cdf of the lower tail
            rhsee_each(z) = rhsee_temp * pdf / ( 1d0 - 2d0*cdf )

        end do
        
        ! compute integeral
        call trapz_integ(z_trapz,rhsee_each,rhsee_derived(t))
                
        cv_ee(t) = rhsee_derived(t)**(-1d0) + psi*hv(t)**(1d0+1d0/nu)/(1d0+1d0/nu)
        error(t) = log10( abs(1d0 - cv(t)/cv_ee(t)) )

        else
        ! Numerical integration when policy is imposed #############################################################
        ! Compute right-hand side of Euler equation using numerical integration ####################################
        ! Then derive Rbar explicitly ##############################################################################
  
        do z=1,n_trapz
                    
            a_next  = rho_a*av(t) + sigma_a*z_trapz(z)
            h_next  = ((1d0-alpha)/psi*exp(a_next)*kv(t)**(alpha))**(1d0/(alpha+1d0/nu))
            y_next  = exp(a_next) * kv(t)**(alpha)*h_next**(1-alpha)
        
            rkhats_next = log ( rbv(t) * (1d0-1d0/lv(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
            rkhat_next  = log ( alpha*((1d0-alpha)/psi)**((1d0-alpha)/(alpha+1d0/nu)) ) &
                        + (1d0+nu)/(alpha*nu+1d0)*a_next - (1d0-alpha)/(alpha*nu+1d0)*log(kv(t))
            xi_next     = exp(rkhat_next)/exp(rkhats_next)
            
            state(1) = log(kv(t))
            state(2) = log(rbv(t))
            state(3) = log(nv(t))
            state(4) = log(xi_next)
            do i=1,nstv
                do j=1,npol
                    x0_next(i,j) = state(i)**expmx(i,j)
                end do
            end do
            x_next  = x0_next(1,:)*x0_next(2,:)*x0_next(3,:)*x0_next(4,:)
            
            if ( xi_next >= 1 ) then
                
                ! Derive l_next and n_next, then check if policy is binding.
                ! If binding, l_next=thre_l and c_next is determined accordingly.
                ! If not binding, use rhsee_next.
                
                pibv_next  = ( exp(rkhat_next) + 1d0-delta )*lv(t)*nv(t) - rbv(t)*(lv(t)-1d0)*nv(t) - nv(t)
                
                if ( pibv_next > 0d0 ) then ! if profit is positive
                    
                    rhsee_next = exp( sum(x_next*eta1_rhsee) )
                    c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
    !                rv_next    = rk_next * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) )
                    
                    n_next = chi1*( chi0*pibv_next + nv(t) ) + n0
                    i_next  = y_next - c_next
                    k_next  = (1d0-delta)*kv(t) + i_next
                    l_next  = k_next/n_next
                    
                    if ( l_next>thre_l .and. n_next>thre_n ) then
                        l_next = (1d0-tau_l)*l_next
                        k_next = l_next*n_next
                        i_next = k_next - (1d0-delta)*kv(t)
                        c_next = y_next - i_next
                        rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) )
                    end if
                
                else ! if profit is negative
                    
                    rhsee_next = exp( sum(x_next*eta3_rhsee) )
                    c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
    !                rv_next    = rk_next * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                    rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) )

                    n_next = chi1*(      pibv_next + nv(t) ) + n0
                    i_next  = y_next - c_next
                    k_next  = (1d0-delta)*kv(t) + i_next
                    l_next  = k_next/n_next
                    
!                    if ( l_next>thre_l .and. n_next>thre_n ) then
!                        l_next = thre_l
!                        k_next = l_next*n_next
!                        i_next = k_next - (1d0-delta)*kv(t)
!                        c_next = y_next - i_next
!                        rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) )
!                    end if

                end if      
                
            else ! if run happens
                
                rhsee_next = exp( sum(x_next*eta2_rhsee) )
                c_next     = rhsee_next**(-1d0) + psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu)
!                rv_next    = ( exp(rkhat_next) + 1d0-delta ) * lv(t)/(lv(t)-1d0) - rbv(t)*lambda
                rv_next    = ( exp(rkhat_next) + 1d0-delta ) * lv(t)/(lv(t)-1d0) / rbv(t) - lambda
                rhsee_temp = beta / ( c_next - psi*h_next**(1d0+1d0/nu)/(1d0+1d0/nu) ) * rv_next
                
            end if

            ! multiply density, and cut off the both tails outside of z_trapz
            call normp( z_trapz(z), dummy1, dummy2, pdf )    ! obtain pdf at z_trapz(z)
            call normp( z_trapz(1), cdf, dummy1, dummy2 )    ! obtain cdf of the lower tail
            rhsee_each(z) = rhsee_temp * pdf / ( 1d0 - 2d0*cdf )

        end do
        
        ! compute integeral
        call trapz_integ(z_trapz,rhsee_each,rbar_derived_temp) ! RHS of Euler, and denominator of Rbar
!        rbar_derived(t) = 1d0/( cv(t) - psi*hv(t)**(1d0+1d0/nu)/(1d0+1d0/nu) ) / rbar_derived_temp
        
        rhsee_derived(t) = rbv(t) * rbar_derived_temp
        cv_ee(t) = rhsee_derived(t)**(-1d0) + psi*hv(t)**(1d0+1d0/nu)/(1d0+1d0/nu)
        error(t) = log10( abs(1d0 - cv(t)/cv_ee(t)) )
        
        end if

        end do
        
        ! End of integration #################################################################################

    error_avg = sum(error(t_drop+1:t_max)) / (t_max-t_drop)
    write(*,*) "EE Error     =",error_avg
    write(*,*) "EE Error max =",maxval(error(t_drop+1:t_max))
    write(*,*) "States when error is max:"
    write(*,*) "k(t-1)=",kv(maxloc(error(t_drop+1:t_max))+t_drop-1),"rbar(t-1)=",rbv(maxloc(error(t_drop+1:t_max))+t_drop-1)
    write(*,*) "n(t-1)=",nv(maxloc(error(t_drop+1:t_max))+t_drop-1),"xi(t)    =",xiv(maxloc(error(t_drop+1:t_max))+t_drop)
    write(*,*) "a(t)  =",av(maxloc(error(t_drop+1:t_max))+t_drop)  ,"pibv(t)  =",pibv(maxloc(error(t_drop+1:t_max))+t_drop)
    write(*,*) "pol_flag(t)    =",pol_flag(maxloc(error(t_drop+1:t_max))+t_drop)
    write(*,*) "max t          =",maxloc(error(t_drop+1:t_max))+t_drop
    write(*,*) "cv(t),cv_ee(t) =",cv(maxloc(error(t_drop+1:t_max))+t_drop),cv_ee(maxloc(error(t_drop+1:t_max))+t_drop)

    open (50, file='from_fortran_to_matlab/ee_error.txt', status='replace', action='write')
    do i=t_drop+1,t_max
        write(50,*) error(i)
    end do
    close(50)

    open (50, file='from_fortran_to_matlab/ee_error_avg.txt', status='replace', action='write')
        write(50,*) error_avg
    close(50)
    
end program
    
! residual function to obtain deposit rate rbv #########################################################
subroutine residual_fn(n_in, temp_in, residual, info_in, t_in)

    use parameter_setting
    use my_toolbox

    implicit none

    integer, intent(in) :: n_in, info_in, t_in
    double precision :: temp_in(n_in)
    double precision, intent(out) :: residual(n_in)
    double precision :: rbv_nls,rks_nls,z_nls,ex_nls,p_nls,pdf_nls,cdf_nls

    rbv_nls = temp_in(1)

    rkhats_next = log ( rbv_nls * (1d0-1d0/lv(t)) * ( 1d0+lambda*(1d0-gamma) ) - (1d0-delta) )
    z_temp     = (rkhats_next - erk_next) / sigma_rk
    z_temp_log = (erk_next + sigma_rk**2d0 - rkhats_next) / sigma_rk
    call normp( z_temp    , cdf    , dummy2, pdf    )    ! obtain cdf and pdf
    call normp( z_temp_log, cdf_log, dummy1, dummy2 )    ! obtain cdf

    ! residual function to solve
    residual(1) = rbv_nls - &
                  ( exp(erk_next + sigma_rk**2d0/2d0) * cdf_log / (1d0-cdf) + (1d0-delta) &
                    - lambda * ( rbv_nls**2d0 * (1d0-gamma) * (1d0+lambda*(1d0-gamma)) * (lv(t)-1d0)/lv(t)**2d0 ) &
                             / ( rbv_nls * (1d0-1d0/lv(t))  * (1d0+lambda*(1d0-gamma)) - 1d0+delta ) /sigma_rk * pdf/(1d0-cdf) )
    
end subroutine